module numeric_func; 
use memory;
use vec;
implicit none;
logical, save :: FUNC_VEC_MOD_TRACE_NEW = .TRUE.;
logical, save :: FUNC_VEC_MOD_TRACE_IF = .TRUE.;
logical, save :: FUNC_VEC_MOD_TRACE_ENTRY = .FALSE.;
!logical, save :: VEC_MOD_TRACE_ENTRY = .TRUE.;
integer, save :: NEW_FUNC_COUNT = 0;


????????
type FUNC_T  !PT - pointer to type; or type-pointer
 real(R14), pointer :: x(:), y(:), dy(:), dy2(:);
 integer mini, maxi;
end type

type FUNC_VEC_T
  real(R14), pointer :: x(:), y(:), drv(:), drv2(:); 
end type
type FUNC_VEC_PT  !PT - pointer to type; or type-pointer
  type(FUNC_VEC_T), pointer :: fp; ! no better way 
end type


contains 
subroutine numeric_func_test();  implicit none; 
  call dbg_info("-->numeric_func_test()"); 
  call func_polynom_test();
  !ok = vec_dot_test_speed();
end subroutine

subroutine func_polynom_test();  implicit none; 
  integer :: nv;
  real(R14) :: minv, maxv, res;
  real(R14), dimension(100) :: px, py;      
  real(R14), dimension(3) :: coeff;      
  call dbg_info("-->new_func_drv_test()"); 

  nv = 3;  minv = D01; maxv = 1.1D0;
!  px => new_step_grid(minv, maxv, nv);   
  call calc_step_grid(px, 1, nv, minv, maxv);
  call dbg_vec("px=", px);
  call assert_max_err("3.px(1)=minv", px(1) - minv, ERR14);
  call assert_max_err("3.px(2)=minv", px(2) - (maxv+minv)/2., ERR14);
  call assert_max_err("3.px(nv)=maxv", px(nv) - maxv, ERR14);
  
  coeff = (/D1, D0, D0/);
  call calc_func_polynom(py, px, coeff);  call dbg_vec("py=", py);
  res = D1;
  call assert_max_diff("1.py(1)", py(1), res, ERR14);
  call assert_max_diff("1.py(2)", py(2), res, ERR14);
  call assert_max_diff("1.py(3)", py(3), res, ERR14);

  coeff = (/D0, D05, D0/);
  call calc_func_polynom(py, px, coeff); call dbg_vec("py=", py);
  res = 0.5D0;
  call assert_max_diff("05.py(1)", py(1), res * px(1), ERR14);
  call assert_max_diff("05.py(2)", py(2), res * px(2), ERR14);
  call assert_max_diff("05.py(3)", py(3), res * px(3), ERR14);

  coeff = (/-D05, -D1, D3/);
  call calc_func_polynom(py, px, coeff); call dbg_vec("py=", py);
  res = -0.5D0;
  call assert_max_diff("513(1)", py(1),res-px(1)+3*px(1)*px(1),ERR14);
  call assert_max_diff("513(2)", py(2),res-px(2)+3*px(2)*px(2),ERR14);
  call assert_max_diff("513(3)", py(3),res-px(3)+3*px(3)*px(3),ERR14);
end subroutine


function new_func2(xp, yp) result(res);  
  implicit none;  
  type(FUNC_VEC_T), pointer :: res;!!! no, do not use FUNC_VEC_PT      
  real(R14), pointer :: xp(:), yp(:);
  integer :: status;
  call dbg_info("-->new_func(xp,yp)", FUNC_VEC_MOD_TRACE_ENTRY);
  allocate(res, STAT=status);
  res%x => xp;
  res%y => yp;
  NULLIFY(res%drv, res%drv2)
  NEW_FUNC_COUNT = NEW_FUNC_COUNT + 1;
  call dbg_int("NEW_FUNC_COUNT=", NEW_FUNC_COUNT, FUNC_VEC_MOD_TRACE_NEW);
  return;
end function
  
  
function new_func(xp) result(res);  
  implicit none;  
  type(FUNC_VEC_T), pointer :: res;!!! no, do not use FUNC_VEC_PT      
  real(R14), pointer :: xp(:);
  integer :: status;
  call dbg_info("-->new_func()", FUNC_VEC_MOD_TRACE_ENTRY);
  allocate(res, STAT=status);
  res%x => xp;
  res%y = new_vec_from(xp);
  NULLIFY(res%drv, res%drv2)
  NEW_FUNC_COUNT = NEW_FUNC_COUNT + 1;
  call dbg_int("NEW_FUNC_COUNT=", NEW_FUNC_COUNT, FUNC_VEC_MOD_TRACE_NEW);
  return;
end function


subroutine del_func(res);  implicit none;  
  type(FUNC_VEC_T), pointer :: res      
  call dbg_info("-->del_func()", FUNC_VEC_MOD_TRACE_ENTRY);
  call dbg_bool("ASSOCIATED(res)=", ASSOCIATED(res), FUNC_VEC_MOD_TRACE_IF);  
  if (ASSOCIATED(res)) then
    call del_vec(res%x);
  end if 
  call del_func_y(res);     
  return;
end subroutine

! keep x
subroutine del_func_y(res);  implicit none;  
  type(FUNC_VEC_T), pointer :: res      
  call dbg_info("-->del_func_y()", FUNC_VEC_MOD_TRACE_ENTRY);
  call dbg_bool("ASSOCIATED(res)=", ASSOCIATED(res), FUNC_VEC_MOD_TRACE_IF);
  if (ASSOCIATED(res)) then
    call del_vec(res%y);
    call del_vec(res%drv);
    call del_vec(res%drv2);
    DEALLOCATE(res);
    NEW_FUNC_COUNT = NEW_FUNC_COUNT - 1;
  end if     
  return;
end subroutine

! note 'y'; meaning only y-values are freed
subroutine set_func_py(oldp, newp);  implicit none;  
  type(FUNC_VEC_T), pointer :: oldp, newp;      
  call dbg_info("-->set_func_py()", FUNC_VEC_MOD_TRACE_ENTRY);
  call dbg_bool("ASSOCIATED(oldp)=", ASSOCIATED(oldp), FUNC_VEC_MOD_TRACE_IF);
  call dbg_bool("ASSOCIATED(newp)=", ASSOCIATED(newp), FUNC_VEC_MOD_TRACE_IF);
  call del_func_y(oldp);  
  oldp => newp;
  return;
end subroutine



function new_func_exp(px) result(py);   implicit none;  
  real(R14), dimension(:), pointer :: px, py;      
  call dbg_info("-->new_func_exp()");
  py => new_vec_from(px);  call dbg_vec("py=", py);
  py = EXP(px);
  return;
end function

function new_func_sin(px) result(py);   implicit none;  
  real(R14), dimension(:), pointer :: px, py;      
  call dbg_info("-->new_func_sin()");
  py => new_vec_from(px);  call dbg_vec("py=", py);
  py = SIN(px);
  return;
end function

function new_func_cos(px) result(py);   implicit none;  
  real(R14), dimension(:), pointer :: px, py;      
  call dbg_info("-->new_func_cos()");
  py => new_vec_from(px);  call dbg_vec("py=", py);
  py = COS(px);
  return;
end function

!function new_func_polynom(px, coeff) result(res);   implicit none;  
!  real(R14), intent(IN) :: coeff(:);      
!  real(R14), pointer :: px(:), py(:);   
!  type(FUNC_VEC_T), pointer :: res;    
!  call dbg_info("-->new_func_polynom()");
!  py => new_vec_polynom(px, coeff);
!  res => new_func2(px, py);
!  return;
!end function

!
SUBROUTINE calc_func_polynom(res, vx, coeff);   implicit none;  
  real(R14), intent(IN)  :: coeff(:), vx(:);      
  real(R14), intent(OUT) :: res(:);      ! vector-x, vector-y
  integer :: i, mini, maxi;
  call dbg_info("-->calc_func_polynom()");
  mini = LBOUND(vx, 1); 
  maxi = UBOUND(vx, 1); 
  call assert_true("LBOUND(res, 1) .LE. mini", LBOUND(res, 1) .LE. mini);
  call assert_true("maxi .LE. UBOUND(res, 1)", maxi .LE. UBOUND(res, 1));
  do i = mini, maxi 
    res(i) = func_polynom(vx(i), coeff);
  end do
end SUBROUTINE


real(R14) function func_polynom(x, coeff) result(y);   implicit none;  
  real(R14), dimension(:), intent(IN) :: coeff;      
  real(R14), intent(IN) :: x;      
  real(R14) :: powx;      
  integer :: i, mini, maxi;
  call dbg_info("-->func_polynom()");
  call assert_true("SIZE(coeff, 1) > 0", SIZE(coeff, 1) > 0);
  mini = LBOUND(coeff, 1); 
  maxi = UBOUND(coeff, 1); 
  powx = 1.;
  y = coeff(mini);
  do i = mini + 1, maxi
    powx = powx * x;
    y = y + powx * coeff(i);
  end do
  return;
end function


end module 